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Abstract 

This is the last in a series of papers on the topological susceptibility in the 
interacting instanton liquid model (IILM). We will derive improved finite tem- 
perature interactions to study the thermodynamic limit of grand canonical 
Monte Carlo simulations in the quenched and unquenched case with light, 
physical quark masses. In particular, we will be interested in chiral sym- 
metry breaking. The paper culminates by giving, for the first time, a well- 
motivated temperature-dependent axion mass. Especially, this work finally 
provides a computation of the axion mass in the low temperature regime, 
m 2/2 = 1.46 10~ 3 A 4 1+ ^° 3 5 °^.48 - It connects smoothly to the high tem- 
perature dilute gas approximation; the latter is improved by including quark 
threshold effects. To compare with earlier studies, we also provide the usual 
power-law m\ = fi^fl)*- ' where A = 400 MeV, n = 6.68 and a = 1.68 10~ 7 . 



1. Introduction 

The purpose of this paper is to continue the study of the interacting instan- 
ton liquid model (IILM) at finite temperature, with an ancillary goal being to 
improve our understanding of the temperature dependence of the axion mass. 
For the first time, we will be able to give a well-motivated axion mass that 
covers all temperatures down to T = 0. 

In [3] we set up the formalism underlying the IILM and developed a numer- 
ical framework to compute the interactions given by an arbitrary background 
ansatz. Simulations at zero temperature were performed with the so-called ratio 
ansatz to determine the parameters that enter the model: the lambda parame- 
ter, A, and the quark masses, m q . 
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Figure 1: The Rjj ansatz does not lead to a well-defined thermodynamic limit at finite temper- 
ature, whereas the Re ansatz does exhibit the correct linear scaling of an extensive quantity. 
Here we display results from quenched simulations, but the same problems persist with dy- 
namical quarks. The instanton number N shows the strongest response to the unphysical 
interactions. 



In this paper we will investigate the IILM at finite temperature based on the 
caloron solution of Harrington and Shepard [14j . Using as input the physical 
parameters we determined at zero temperature, we will study chiral symmetry 
restoration, based on the ideas of instanton-anti-instanton molecule formation 
151 ] . 16j . [3a |. and determine the topological susceptibility. 

In order to deal with light, physical quark masses, we study the thermo- 
dynamic limit. As mentioned in the previous paper ; 
interactions derived previously in [3^ are deficient in this respect. 



we found that the 
Although 

the old and new interactions agree rather well at zero temperature, this is no 
longer the case at T ^ because of an unphysical behaviour for the instanton- 
instanton interaction that onl y d ecays as 0(l/R s ), with R s the spatial instanton 
separation, see Eq. (3.11) in [39(. This long-range interaction prohibits a ther- 
modynamic limit as it is not integrable, see Fig. [TJ In their paper, the authors 
do discuss this long-range interaction and report they found the 0(1/ R s ) dyon- 
dyon behaviour for a wide range of intermediate separations. It might well be 
that for the simulation boxes used in a subsequent numerical investigation (see 
34j) the interactions are still well described by this ansatz, i.e. that the spatial 
extent of the box is bounded by these intermediate separations. However, for 
studying the large volume behaviour this ansatz is not appropriate. Removing 
this particular part of the interactions, we were able to retrieve a well-defined 
thermodynamic limit. 

Apart from this deficiency, it seems obvious that at finite temperature it will 
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be much harder to find a good parametrisation for the action of the background 
ansatz because the underlying 0(4) symmetry is broken and the constituent 
gauge fields are more complicated. This fact was already pointed out in_|39(. 
We improve the existing interactions by extending the formalism set up in [44] to 
finite temperature; as we will see, the only difficulties are of a technical nature. 

With this paper we will achieve our initial aim of computing the temperature 
dependent axion mass. So far the mass was computed within a dilute gas ap- 
proximation which breaks down at low temperatures. The connection between 
the high temperature regime, where the dilute gas becomes ever more accurate, 
and the zero temperature result, obtained through chiral perturbation theory, 
has been performed in a rather crude manner up to date: either by unsmoothcd 
matching 41 , i[ or by an ad hoc interpolation prescription [22[ • Our determina- 
tion of the axion mass will for the first time give a well-motivated interpolation 
between the zero and finite temperature regimes. Comparison with lattice data 
will allow for a critical evaluation of the systematic uncertainties. In partic- 
ular, considerations regarding the anthropic axion with large decay constant 



27. 42, 40] are potentially very sensitive to the non-perturbative effects of the 
QCD phase transition, when their mass becomes sizable. 

In section [2] we will re-derive the finite temperature interactions for the ra- 
tio ansatz. We will then discuss, in section [31 the new elements that finite 
temperature introduces in the numerical framework. After these technical pre- 
liminaries we will have a short investigation of the topological susceptibility in 
the quenched sector in section |4l before we discuss the main numerical results 
of the unquenched IILM regarding chiral symmetry restoration, section [3 and 
the topological susceptibility and axion mass, section [6] 



2. Interactions in the IILM at finite temperature 



In terms of the 't Hooft potential l+II, the Harrington-Shepard caloron [14j . 
an infinite sequence of singular gauge BPST instantons [f| along the Euclidean 
time direction, is given by 



A a = _ n ab r b d v IL(x,{y, P }) 



where O is the colour matrix in the adjoint representation of the embedding 
SU{2) — > SU(3), and = fj b = rf ) for instantons (anti- instantons) ; r\ 

are the 't Hooft symbols. The 't Hooft potential has the following form 

2 S inh 2gE 

U(x,{y,p}) = ^. g 2gt , (2) 

pr cosh — cos 

with r 2 = (x — y) 2 and t — x 4 — y 4 ; the collective coordinates are: y the centre, p 
the size and O the colour orientation. At finite temperature bosonic quantities 
such as n are periodic in the Euclidean time direction with period /3 = l/T. 
Note that n approaches the zero temperature instanton potential in the singular 
gauge for f} — > oo. 
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We will use as background the ratio ansatz and, as for the zero temperature 
case, only consider two-body interactions. The gauge field is then given by 

with O = 0[ O2. This pair interaction has been derived in [3^|, and we will 
refer to it as Rh; the corresponding forces derived in this paper will be denoted 
by Re- However, in studying the volume dependence of various quantities, 
we noticed that these interactions did not allow for a thermodynamic limit. 
The problem can be traced back to the log term in the instanton-instanton 
interaction, (3.11) in [39}, only decaying like 0(1/ R s ) for large R s , where R s 
is the spatial separation of the pair. Note that this term is attributed to the 
dyon-dyon interaction for intermediate separations, /3 <C R s <C p 2 //3, in the 
high temperature limit. 

The ratio ansatz has the same functional form in terms of II as for zero 



temperature, and so we can use our previous result (44| to write 



F^Kv = I+ ( Tr0t ° + (f)Ov)^u)J + (f)Or,)p» P uI»» 

The different terms are given in appendix |Appendix A| Due to charge renor- 
malisation the action, S[A] — j% J F^ V F^ V , acquires a quantum contribution 
in the form of the running coupling constant. The classical interaction is given 
by 

SyS = V 12 = {S[A]/S -2), (5) 

where So = 8Tr/g 2 is the single instanton action. The quantum effects substitute 
g in So for the running coupling constant; the RG scale is estimated by the 



geometric mean y/p~ip2, as proposed in [38|, |34 1 . 



First, we will look at the dependence of the interactions on R t , the instanton 
separation in the (imaginary) time direction. Compared to the zero temperature 
case, the differences are substantial even for pairs with equal sizes. The major 
difference comes from the fact that the Rh interactions are not periodic. For 
unequal sizes, the difference is even more pronounced as was the case at zero 
temperature, see Fig. [2] The reason is again that the functional form on the 
instanton sizes in Rh is not general enough. 

The interaction between instantons and anti-instantons does not lead to as 
close a match between Re and Rh as in the T = case, but is still qualita- 
tively similar for not too large separations. In the instanton-instanton case, 
however, we have found significant differences, see Fig. |3] we do find the dyon- 
dyon behaviour for intermediate distances but for separations R s p 2 //3 the 
functional dependence changes into an integrable 0(Rj 4 ) for Re whereas the 
non-integrable dyon-dyon interaction persists in the Rh ansatz. That the large 
separation fall-off should be different from the dyonic regime is clear from the 
discussion in [lij where the authors show that the caloron field develops a dipole- 
like character in the far-field region. Note that this leads to a three-dimensional 
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Figure 2: The major difference follows from the fact that Rh interactions are not periodic in 
Rt. This is clearly a deficiency of the analytical formulas because it follows directly from (J4j) 
that the interactions should have period j3. For spatial separations, the main differences occur 
for unequal size parameters, e.g. p\/p2 = 3 in this case. The reason is that the dependence 
on the sizes is more complicated than the functional form, ^Jpxpi, used in Rh- (We have set 

P = \IpI+pI-) 
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Figure 3: The instanton— instanton interaction behaves very differently for Ru and Re- In 
the high temperature, large size limit we do see the dyon— dyon behaviour but for larger 
separations, R s S> p 2 /P, the interaction decays much faster, 0(RJ 4 ). This fall-off behaviour 
is integrable and we can study the thermodynamic limit, in contrast to the purely dyonic 



interaction. (We have set p = 




dipole-dipole interaction for an instanton-anti-instanton pair, but that the Rh 
ansatz retains the zero temperature (four-dimensional) dipole-dipole intcrac- 
tiorJE 

The fermionic interaction follows from the quark wave function overlaps 

{p + m) tf = (& \p + m%) = pij + mS l3 . (6) 

Even though the set of eigenf unctions {^} is generally not an orthonormal 
basis, the extra contributions are neglected; effectively, we treat as being 
orthonormal which leads to the diag onal mass term. The fermionic zero mode 
at finite temperature is given by 13 1 , 12 1 , [39[ , 

cos 2* 

* = "co^' (9) 



1 instanton— instanton pairs decay slightly faster, at 0(1/ R^), however, and their interaction 
is not of dipole type; note that this happens at T = too, where instanton-anti-instanton 
pairs exhibit a dipole-dipole interaction but instanton-instanton pairs don't [jj. 



G 



with ip aa — e aa , normalised according to £12 = 1, and Ui the collective coordi- 
nates for the colour embedding in the fundamental representation. 

The overlaps Tja = J ^jip^A have a slightly more complicated form than 
their T = counterparts and are given by 

Tia = I — ^ (\Tr{UT+)Ip - l -Tz{UT+T a )% a J^ a 

jRZxS 1 ^ Pi PA V 2 2 

+ l -Tr(Ur a T+)^ a K^ }j . (10) 



The different contributions can be found in appendix |Appcndix B 



The difference between the Re and Rh interactions is quite large for tem- 
poral and spatial separations, see Fig. |4] However, it is not straightforward 
to compare both ansatze because the Rh overlaps are computed on a different 
background [39| . the sum ansatz, whereas we use the full ratio ansatz in Re- 
As in the gluonic case, the quark overlaps are not periodic for the Rh ansatz 
as opposed to those of Re- As always, unequal sizes increase the differences 
between the Rh and Re ansatze even more. For large separations, when the 
ratio ansatz becomes indistinguishable from the sum ansatz, and for equal sizes, 
we find very good agreement between Rh and Re- 

The total interaction, after normalising to the dilute gas approximation, is 
given by 



pairs 



c , x , lndet(I+^) ,Q<0 

so(vPiPjm 3 -^\ lndet(I+ ^ ) )Q>0 , (id 



where Q = Ni — Na is the topological charge and Nf the number of quark 
flavours. For details see 

Even though the analytic formulas of the Rh ansatz are not periodic, this 
is not really a major shortcoming because periodicity can be realised by folding 
back the instantons into the fundamental interaction region R t € [— /3/2, /3/2]. 
We demonstrate this for the total interaction in Fig. \5\ We can clearly see that 
the Re ansatz is intrinsically periodic. 



3. Numerical Implementation 

3.1. Interpolation and asymptotic matching 

As for zero temperature, we will need to set up a grid for the numerical 
evaluation of the two-body interactions. The look-up tables will depend on four 
variables: the spatial separation R s , the temporal separation R t and the two 
sizes pi and p2- The colour degrees of freedom O, or equivalently U, have been 
completely factored out and can be treated exactly. 

As we have seen in section [1] the dyon-dyon regime is characterised by a 
fairly slow fall-off. However, not all bosonic interactions in appendix |Appcndix A.l| 
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Figure 4: The large discrepancy between Re and Rh for both temporal and spatial separations 
is due to the fact that Re uses the full ratio ansatz in the Dirac operator whereas for Rh a 
sum ansatz is used. Periodicity, which follows directly from {((J, is lacking for overlaps in the 

Rh ansatz but is realised in the Re ansatz. (We have set p = W p\ + p\.) 
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Figure 5: Folding back the temporal separation into the fundamental interaction region, 
Rt £ [— /3/2, /3/2], we explicitly retrieve periodic interactions for Rjj- The Re ansatz is 
intrinsically periodic. 



decay this slowly; rather they can be grouped according to polynomial and ex- 
ponential decay. The fermionic overlaps fall-off exponentially, but slightly more 
slowly than their bosonic counterparts. Therefore, we found it advantageous to 
define three different grids: 

• For the polynomial grid, the maximal separation i?" lax will depend on 
p = irp 2 /(3, the natural size parameter for separations beyond (3. Note 
that for very small sizes i?" lax < (3 the prescription for i?™ ax should 
switch over to the T = case; this is implemented by setting i?™ ax = 
max(a!g/9, a^irp 2 / j3). 

• The exponential decay sets in at R s ss /3. To accommodate very small 
instantons we set again i?™ ax = min(aQ p, afp /3). 

• The grid for quark overlaps is set in an analogous fashion, i.e. i?™ ax = 
min(ag p, /3). 

The constants a? are fine-tuned so as to achieve fast and stable numerical in- 
tegrations and a good matching to the analytic expressions used for separa- 
tions beyond J?™ ax . In the temporal direction the grid is bounded by |i? t | = 
f3/2; care needs to be taken again for small instantons and we set |-R™ ax | = 
min( J R™ ax ,/3/2). 

The size distribution is supported on 30 grid points and for each pair (p.;, pj) 
the R s —Rt plane consists of 30-29 nodes; this leads to a total of roughly 1 million 
interpolation points. The grid is 'logarithmic' in the size and R s direction. A 
typical R s — R t plane is shown in Fig. [5] 
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Interpolation among four points 




Figure 6: It turns out that the interactions are more or less identical on the R s — Rt plane, 
when expressed in units of p\ + p|. We therefore choose the grid to be irregular, i.e. the 
points on the R s — Rt plane are different for different values of the sizes. Since instantons 
are localised field configurations, the separation between grid points grows geometrically in 
the R s direction. In the p directions the mesh is regular, but still 'logarithmic'. The inner 
and outer regions are matched to asymptotic expansions: for small separations the matching 
is performed radially from the origin because in this region the interactions are roughly 0(4) 
symmetric; for large separations we match along the R s direction with constant Rt because 
the interactions become t-independent. (Note that for graphical reasons not all 870 grid 
points arc displayed and that we actually display log R s . Also the actual grid only contains 
the R a > points.) 
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Figure 7: In order to derive simple analytical expressions in the limit of small and large pair 
separations, we approximate the 't Hooft potential to have a simple functional form in different 
integration regions. This procedure works well for instantons with sizes that do not exceed /3. 
As we will see the plasma screening effects limit instanton sizes to be rather small compared 
to 0. 



The minimal and maximal sizes supported by the grid are fixed as in the 
T = case, i.e. they are chosen so small, respectively large, that they correspond 
to a very low, respectively high, quantile of the normalised instanton density. As 
we will see in section 13.21 at finite temperature the instanton density becomes 
temperature dependent, reflecting the screening of coherent background fields. 
Instead of computing a grid for every temperature we run simulations at, we 
can exploit the following scaling transformations: under coordinate rescaling, 

X ^ 7 OCX 

V( Pi ,Rj,p) = V(jH/P,Ri/0), (12) 
T{fH,Rj,P) = h?(jHlP,Rj/P). (13) 

These can be used to transform the grid and interactions to different temper- 
atures. We just need to make sure to choose the original sizes large enough to 
accommodate the low temperature behaviour, i.e. p max = A. We defined the 
grid and the interactions at /3 = 1 with p e [0.01, 1.6]. 

As for the zero temperature case, the matching consists in deriving asymp- 
totic formulas / asy that are patched on to the numerically integrated interactions 
according to 

m = fasy (R)M^}- (14) 



where i? m is the matching point. For details see [44|. At T =^ 0, however, we 
match differently for small and large separations. In the former case the inter- 
actions will be approximately 0(4) symmetric and the matching is performed in 
an 0(2) symmetric way on the R s — R t plane, i.e. along a ray connecting the ori- 
gin with (R s ,R t ). For large separations the interactions become t- independent 
and we match along a line of constant R t . 

To derive the large separation asymptotic formulas, we consider two different 
integration regions, see Fig. [7J Therefore, we effectively approximate integrals 
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Figure 8: The instantons Ii and I2 are so far apart that within the shaded region that give 
the dominant contribution to the field strength of each, the other's field strength is roughly 
constant and fixed at x M — R^ —Ru,. We can then safely extend the integration region to 
be all of S 1 X R 3 , with a negligible error due to the rather strong localisation of the individual 
instantons. The functional form of the relevant 't Hooft potential changes from II = p 2 /r^ d 
to n = np 2 /l3r 3d at r 4d = 0/2, see Fig.[3 



by 

1*1 +P I , (15) 

JSixR 3 J.S 3 x [0./3/2] J S 2 x [,3/2,oo] 

where we exploit the respective spherical symmetry of the 't Hooft potentials. 
In both regions we only have to deal with rational functions which can be easily 
integrated exactly. Apart from this extra complication the strategy is the same 
as in the zero temperature case and is illustrated in Fig. [BJ The results are given 
in appendix I Appendix A.2T| and |Appcndix B.2.1| 

The main objective for small separations is to capture adequately the sin- 
gularity structure, for which it is sufficient to restrict ourselves to the T = 
region. To capture contributions that do not blow up, we add the large sepa- 
ration contributions truncated at the pair separation (i? t ,i? s ). The results are 
given in appendix | Appendix A.2!2| and |Appcndix B.2.2| Apart from considering 
different functional forms for the 't Hooft potential this is the same procedure 
as for the T = asymptotic interactions and is illustrated in Fig. 

As in the T = case, the quark overlaps in the small separation region are 
qualitatively incorrect. However, this region is again dominated by the gluonic 
repulsion. The latter is not as well approximated as in the T = case but the 
agreement is still on the 10% level and, most importantly, qualitatively correct, 
see Fig. 1101 In contrast, the large separation asymptotic formulas work well. 

3.2. Biased Monte Carlo 

The IILM is defined by its partition function, 

oo Ni N A 

Nj,N A A ' i j 

with S'int defined in (fTTj) and the single instanton density given by d(p) = 



do(p)dr(p)- The zero temperature contribution do is given in [44j. The fi- 



nite temperature term describes the screening in the plasma of coherent field 
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Proper small separation asymptotic integration down to R 




Large separation asymptotic integration up to =| 



Figure 9: The instantons I\ and I2 are strongly overlapping. We approximate the integral 
by first integrating over Ii keeping I2 fixed at R^ as in the large separation case, see Fig. [8] 
but with upper limit R/2; to this we add the analogous contribution from 12- The possibly 
singular behaviour is picked up by integrating from infinity down to R and approximating the 
arguments by x M — R^/2 PS x M and x M + -R M /2 !=s x M , respectively. 



excitations exceeding the inverse temperature scale /3, and is given by [12[ [34 1 



dr(p) = exp 



--{2N c + N f ){vpT) 2 



N c N f 

1 + — i- 

6 6 



- ln(l + (tt P T) 2 /3) + 



0.15 



(1 + 0.15(7rpT)- 3 / 2 ) 8 



(17) 



Apart from book-keeping related technicalities, the major challenge is to ad- 
equately simulate instanton-anti-instanton molecules, the structures responsible 
for the chiral symmetry restoration within the IILM [lj|, [ll|, [35l. Detailed 
numerical studies at finite temperature have been performed in [34j j . For the 
mass parameters used in these latter simulations no technical problems were 
encountered. However, for the small quark masses that we determined in the 
T = simulations, we found that standard Monte Carlo techniques faced severe 
problems with the strong attraction between instanton-anti-instanton pairs. In 
[45| we have argued that this technical problem will occur quite generally for 
the semi-classical expansion with small quark masses. In that paper we have set 
up the framework to deal with these algorithmic issues by adapting techniques 
from chemical engineering and computational chemistry and developed for the 
study of strongly associating fluids. 

The instanton-anti-instanton pair interactions become ever stronger and 
more localised for small quark masses, and random sampling methods will gener- 
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Figure 10: The small separation asymptotic approximation is only qualitatively correct. It 
overestimates the correct interaction for moderately small separations because it underesti- 
mates the quark contribution. For very small separations, where the exact quark interaction 
is indeed negligible, the gluonic approximation underestimates the exact result. However, the 
grid covers rather small separations so that the asymptotic interactions for strongly overlap- 
ping instanton— anti-instanton pairs are rarely needed. In those cases where it is needed wc 
get a qualitatively correct repulsion. (The temporal separation Rt has been chosen very small 
in order to see the repulsion). 
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ically miss these configurations. We therefore need to preferentially sample the 
attraction centres. This can be achieved by resorting to biased Monte Carlo 
schemes. They exploit the large redundancy in devising transition probabilities 
Pij that satisfy detailed balance, 

prPi^pppji, as) 

where the equilibrium distribution P eq is given by the partition function. Specif- 
ically, the degeneracy follows from the fact that the transition probability con- 
sists of two parts, the proposal probability, Vij, and the acceptance probability, 

Pij = VijAij . (19) 



To satisfy (|18|) . the acceptance probability is given by 



Aij — min 



(20) 



which is the Metropolis prescription. 

Recently, efficient and general purpose algorithms have been developed to 
achieve this importance sampling; we will use 48], the Unbonding-Bonding al- 
gorithm (UB). This algorithm focuses not on the union of all the interaction 
regions, a complicated and case-specific geometrical problem, but on the indi- 
vidual interaction regions and all the possible routes that lead to the same final 
state. 

For the IILM, we choose the bonding region to consist of the two discon- 
nected valleys of strong interaction, see Fig. [11] and take them to be hypercubes 
for simplicity. The precise size of the box is a free parameter and will be chosen 
to achieve the fastest convergence. We will check that the results do not depend 
on the bonding box parameters, within statistical uncertainties. We will now 
present the transition probabilities for the different types of updates; details can 
be found in jljj . 

The forward and backward transition probabilities for the UB moves of in- 
stanton i are given by 

V *= E p&jj+W'. ( 21 ) 

3 

Nf{i) 
3 

(23) 

with Sf = 1 if i is bonded and Sf = otherwise. We denote by Nf (N^) the 
number of bonded instantons (anti-instantons) and Nf(i) is the number of anti- 
instantons that instanton i is bonded to, and analogously for N^(i). Unprimcd 
quantities are evaluated before the move whereas primed ones denote the same 
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Gauge singularity/repulsion Bonding Box 




Strong and short-ranged R s / p 
interaction 



Figure 11: The strongest interaction is achieved for the colour orientation U = I and located 
in the two troughs. For simplicity we will choose the bonding box to consist of the two 
hypercubcs displayed by dashed lines. 



quantity after the move. The individual bonding and unbonding transition 
probabilities are given by 



P &'*)- NjNaVj' (24) 

V * = WV- (25) 

The bonding move consists of choosing uniformly an instanton and an anti- 
instanton, and to place the instanton i uniformly in the bonding region Vj of 
the anti- instanton j. The unbonding move consists of choosing uniformly one 
of the bonded instantons and to place it uniformly in the simulation box. 

Insertion and deletion will be constructed along the lines of the UB algorithm 
by either placing the instanton i into the bonding region of an anti-instanton or 
removing the bonded instanton i. Adding up all possible paths, including the 
unbiased one, we get 

K B d) 1 x 1 

V NlN -= Pb £ +(1- Pb )-, (26) 

3 J 

VN- I N 1 =Pb§r § + ^-Pb)^T i - (27) 
Here pb is the a-priori-probability to perform biased moves. 
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It has also been argued that cluster moves are important to achieve good 
mixing of the Markov chain [28] . We assume that I — A pairs are the dominant 
clusters that form. The pair displacements consist of translating the pair as 
a whole or to displace one of the constituents within the bonding region of 
its partner, and the resulting transition probabilities are obvious. The pair 
insertions and deletions are given by 

Ci„ 1/1 i \ ill 

7V ~ = siff + 2 w, w A ■ (29) 

where Np is the number of pairs. The last terms in each line follow from ran- 
domly and independently inserting/deleting an instanton and an anti-instanton. 
We found that inclusion of these unbiased moves enhances the acceptance rates 
at high temperature when we approach the dilute gas limit. 

So far we have only discussed the spatial arrangements of the instantons. 
However, the interaction also depends on the colour degrees of freedom (To| . 
This further technicality can be dealt with by adapting techniques for orientation 
dependent forces in molecular dynamics simulations. To this end we define 
the measure in colour space by the Boltzmann factor of the pair interaction. 
This ensures that U is chosen so as to increase the interaction. More precisely, 
given the position of the instanton within the bonding box of its partner anti- 
instanton, we want to sample the following equilibrium distribution: 

exp ( -Sf^iU!, Ua,xi,xa)) 
V{U I \x I ,x A ,U A ) = ±— (30) 

/ dCZ/cxp (-Sf^"(Ul,UA,Xl,X A )j 

where dUi is the Haar measure over SU(3). Note that we neglect the influence 
of the neighbouring instantons and anti-instantons; we could include them but 
since the fermionic interaction involves a determinant the computation would 
become rather costly. 

We cannot sample this distribution analytically, but we can sample it exactly 
within a Monte Carlo scheme [lOj. Using the fact that the Haar measure is 
invariant under group composition and also that the interactions (UJ) and (JTUJ) 
depend on U = U\Ua-, it is natural to work with U. To sample ([5T)]). we choose 
a set of N different colour matrices uniformly over SU(3), i.e. according to 
the Haar measure. We select one among those N matrices according to the 
probability 

e-sr 1 

Pi = ZTiv S5=iT > ( 31 ) 



e 



with Sf 3 "" = Sf nt (Ui, X/, xa)- Now, remember that the transition probability is 
given by a proposal and acceptance probability, and that the former includes any 
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a-priori-probabilities for generating the proposed state. Thus, the orientational 
bias includes the generation probability of the set of matrices, and is given by 



Vi = T , N Pl . (32) 

The selection probability and the first N volume factors are obvious. The ad- 
ditional N — 1 volume factors need further explanation. Given the form (|20[) 
for the acceptance probability, apart from the N matrices used in (|31[) , we need 
to generate another N — 1 trial matrices that combine with the matrix selected 
in the forward move to form the set of N matrices necessary to evaluate the 
orientational bias, (|3ip . in the backward move. Finally, the extra enhancement 
of N is necessary because we are not so much interested in Ui being generated at 
the ith trial but rather that the value Ui is selected. Thus we must marginalise 
over the label i: adding up all the permutations of the set of trial matrices mod- 
ulo the permutations of the subset of the N — 1 matrices that are not selected 
produces the extra N factor. Of the 2N — 1 volume factors 2N — 2 will cancel 
in genera]! and so we can simplify (|32|) to 

e -Hi 
N l^j e 

which we recognise as an approximation to (|30[) : note, however, that within the 
MC scheme is exact, i.e. the results do not suffer from any discretisation errors. 
Including the orientation bias amounts to multiplying the bonding box volume 
factors Vj by (1531 . 

Wc found that the interaction are rather self-similar when expressed in units 
of p = \J pj + p\. Thus, the bonding region depends on both particles that 
make up the pair through a function of their sizes. We choose the functional 
form of the size dependence to be given by 

V ] ^V lJ ^V VB (p 2 l +pf j 2 . (34) 

By construction, V\jb is independent of the specific pair, and we will tune it to 
achieve fast convergence. 

The IILM does not lead to confinement; from a heuristic point of view this 
follows from the simple fact that the Harrington-Shepard calorons have a van- 
ishing Polyakov loop, the order parameter of the phase transition. It has been 
in argued in [37} that at low temperatures the instanton density should only 
depend very weakly on temperature and that, in particular, the result (JTTJ) is 
only applicable in the plasma, i.e. deconfined, phase. In the hadronic phase 
the fundamental degrees of freedom are not screened and the heat bath consists 



2 One might think that all 2N — 1 volume factors drop out, but in general we might add 
new moves that select U with a rule different from (13 1 H ; in such a case we only need 2N — 2 
trial matrices. 
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Figure 12: The screening factors l|35|l lead to constant mean instanton sizes in the confined 
phase, below T„, as has been found from lattice investigations [jj. As expected, the factor d.2 
leads to a slower and more symmetrical switch-on of the screening. 

mainly of pions. We will follow the prescription put forward in [34j and include 
in (fTT|) a phenomenological term that mimics the transition from the confined 
to the deconfined phase. We have decided to investigate this with the following 
two choices, see Fig. [12] 

d 1 T (p)^d (p)d T (p)^ 7 

d 2 T (p) = d (p)(l + a T (dT(p)-l)) , 

1 / T — T \ 

a T = -(l + tanh^^J . (35) 

The functional form d\, for the transition between hadronic and plasma phase 
has been used in [34j. The functional form d\ is chosen on the basis that, 
intuitively, it should lead to a more symmetrical behaviour around T». With 
these two choices we will get a rough idea on the systematics due to these 
phenomenological terms. We've introduced two new parameters into the IILM; 
they will be fixed by comparing our results to available lattice data. 

So far these modified screening factors are purel y p henomenological. How- 
ever, the more general non-trivial holonomy calorons [23, HH 0| [Ifij], that might 
play a crucial role in driving the confinement/deconfinement phase transition, 
have been shown to lead to mean instanton sizes that depend very weakly on 
temperature in the confined phase [111 ]. Given the assumptions of that paper, 
the change of the Polyakov loop with temperature is responsible for the differ- 
ent screening behaviour of instantons in the confined and deconfined phase: in 
the former, the equilibrium state corresponds to an ensemble of maximally non- 
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trivial holonomy calorons which are not screened, whereas the latter consists of 
trivial holonomy Harrington-Shepard calorons which are strongly screened. We 
can therefore interpret the modified screening factors (|35l) as an effective de- 
scription of the change of the Polyakov loop as the system evolves through the 
phase transition. Generalising the IILM to include these calorons in the future, 
we might be able to describe the plasma screening effects self-consistently. 

Note that for low temperatures the IILM should not be very different from 
the T = case, which was well dealt with by ordinary Monte Carlo. We 
can use this regime to test the more advanced Monte Carlo scheme against 
unbiased simulations. In Fig. [T3] we plot the thermalisation history at such 
low temperatures for two different set of quark masses; the latter are given 



by the physical masses determined in [44( and the second set has the same 



dimensionless strange quark mass and identical light quarks with a value of ten 
times the physical up quark. We can clearly see that the histories for the biased 
and unbiased simulations are very similar. 

As the temperature is increased, ordinary Monte Carlo fails for the set of 
physical quark masses. For intermediate temperatures, around T*, there is a 
substantial difference in the results from biased and unbiased simulations. Even 
for the set of larger masses there are differences visible to the naked eye, see 
Fig. 1141 For these temperatures instanton-anti-instanton pair formation, the 
mechanism that drives the chiral symmetry restoration in the IILM, is supposed 
to be very important. We find that the number of pairs is indeed quite different 
with (Np) « 20 for biased and (Np) pa 7 for unbiased simulations and physical 
quarks. For the larger quark masses we find that the difference in the number 
of pairs between biased, (Np) pa 16, and unbiased, (Np) pa 11, simulations is 
less pronounced, as expected. 

At temperatures well inside the plasma phase, ordinary Monte Carlo fails 
completely to describe the system accurately for physical quark masses, see 
Fig. 1151 In this case the number of pairs is very different, (iV£ iased ) Pa 18 and 
(N p nhiased ) w 1. Note, however, that the number of pairs has decreased for the 
higher temperatures. 

Simulations at larger masses converge (slowly) to the correct equilibrium dis- 
tribution. Actually, ordinary Monte Carlo perform better at high temperature 
than around T* for these larger quark masses. The reason is that the ensemble 
equilibrates in a very dilute state, so that there are only a few instanton-anti- 
instanton pairs, (7V£ iasod ) pa 4.5 and (N P nh ' msod ) pa 1.5. Note, again, that the 
pair concentration has dropped with respect to temperatures around T* . It turns 
out that, quite generally, the molecule concentration drops with increasing tem- 
perature, because the energy-dominated pair configurations have a smaller and 
smaller entropy and are outweighed by the entropy-dominated configurations; 
the latter correspond to a truly random ensemble of instantons at a high enough 
temperature. This will happen for physical quark masses as well, albeit at higher 
temperatures. 

To be a bit more precise, consider the contribution to the partition function 
of an instanton-anti-instanton pair, which we approximate by its dilute gas limit 
and by the contribution from the bonding box. The latter dominates, and hence 
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Figure 13: In the hadronic phase, i.e. without the plasma screening factor l|17|l . biased and 
unbiased simulations agree very well. This is to be expected because we know from the 
T = simulations, see |44| . that ordinary Monte Carlo leads to good sampling. These low-T 
simulations serve thus to check the biased Monte Carlo scheme. 
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Figure 14: Around the confinement /deconfinement transition, modeled by H35H , the formation 
of instanton-anti-instanton pairs is important as it drives the chiral symmetry restoration in 
the IILM. We can clearly see that for the small, physical quark masses (top) ordinary Monte 
Carlo fails to sample the ensemble correctly. For the large quark masses the effect is much 
less pronounced but there is a clear systematic difference in the mean instanton number. It 
can be attributed again to fewer pairs. 
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Figure 15: In the plasma phase, ordinary Monte Carlo fails drastically for the physical quark 
masses (top). This trend will, however, not persist to higher temperatures. Depending on the 
quark masses the energy-dominated pair configurations will be outweighed in the partition 
function by entropy-dominated configurations at sufficiently high temperatures. For larger 
masses this happens at lower temperatures (bottom). 
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biased simulations might be needed, if 



1 < 



d{p) AUexp (-S fat ) 



(36) 



The contribution from the two instanton measures have cancelled on both sides, 
and we used V = N/2 J d for the dilute gas. The overbar is meant to indicate 
the average over the bonding box, size and colour orientations. 

Using the 1-loop formula for the instanton measure and choosing one specific 
size with Tp = const, this leads to 



We only display explicitly the temperature dependence for the one-loop in- 
stanton measure (first term), the gluonic interaction (second term), the quark 
interaction (third term) and the scaling with temperature of the bonding box 
(fourth term). It turns out that the gluonic contribution is bounded by 5 for 
Nf — 3. Hence the logT term is negative and the dilute gas approximation 
becomes better and better as the temperature increases, see also Fig. [TH1 For 
higher temperatures, when more quark flavours become active, the trend seems 
to be reversed. Note, however, that this conclusion completely neglects the 
high-frequency quantum interaction, which should be investigated in the future 
to study this issue further. 

The UB algorithm is only useful if the results do not depend strongly on the 
precise implementation, i.e. the actual definition of the bonding box. We checked 
this by using different bonding boxes for three different temperatures, i.e. above, 
around and below T* , see Table [TJ We found that the results agree very well for 
temperatures below T*. This corroborates our expectations that biased Monte 
Carlo is not essential in that regime. For a modest sample of 200 independent 
configurations^), we also find rather good agreement for the higher temperatures 
where random sampling fails. We have also checked the dependence on the a- 
priori-probabilitieg_j we find again that the results only depend weakly on these 
choices for the small sample size we have used. Autocorrelation times, however, 
depend much more strongly on the bonding box; this allows us to fine-tune the 
parameters to achieve efficient sampling. The bonding box B4 will be used for 
the final simulation. 

3.3. Fermionic determinant 

We argued that for sufficiently high temperatures, the IILM goes over to a 
dilute ensemble. It turns out that for high dilution the numerical manipulation 



3 This relates to (N); the autocorrelation time for (Q 2 ) was always smaller, effectively 
leading to a larger sample size. The better agreement on the latter quantity corroborates the 
expectation that the differences among bonding boxes vanish for infinite sample size. 

4 The data in Table [T] corresponds to one particular choice for the a-priori-probabilities. 
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Figure 16: We display a (spatial) snap-shot of a typical configuration. For low temperature 
(top) spatial correlations are not very pronounced and the system equilibrates in a 'random' 
state. In the region T ~ T* (middle) a higher concentration of instanton-anti-instanton pairs 
can be seen, signalling the restoration of chiral symmetry in the IILM. At higher temperatures 
(bottom) the ensemble becomes dilute: the energy gain in pair-formation is outweighed by 
the large entropy gain of a random distribution of instantons throughout the box. 
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Bi 


B 2 


B 3 


B 4 


£5 


B 6 


(N) 


99(1) 


100(1) 


97(1) 


102(1) 


98(1) 


98(1) 


(Q' 2 ) 


0.81(1) 


0.81(1) 


0.81(1) 


0.78(1) 


0.81(1) 


0.779(9) 


(uu) 


230(9) 


219(8) 


222(9) 


217(9) 


211(9) 


225(8) 


(N) 


105(1) 


101(1) 


96.7(9) 


99.3(8) 


98.6(9) 


95.3(9) 


(Q' 2 ) 


0.36(1) 


0.31(2) 


0.39(2) 


0.30(2) 


0.347(8) 


0.349(6) 


(uu) 


84(8) 


68(6) 


80(8) 


69(7) 


75(7) 


72(7) 


(N) 


102.0(9) 


106(1) 


99.2(9) 


106(1) 


108(1) 


102.2(9) 


(Q 2 ) 


6.9(3) 


7.0(3) 


7.0(2) 


6.7(3) 


6.8(2) 


6.9(1) 


(uu) 


1320(30) 


1300(30) 


1300(40) 


1200(30) 


1310(40) 


1210(30) 





Bi 


B 2 


B 3 


B 4 


B 5 


B 6 




0.4 


0.4 


0.4 


0.6 


0.23 


0.17 


Jbb 
i,min 


0.2 


0.25 


0.3 


0.2 


0.35 


0.4 


Jbb 


2 


1.8 


1.6 


2 


1.2 


1 



Table 1: The temperature increases from top to bottom. We find rather good agreement 
between the different bonding boxes. The sample has a modest size of 200 independent 
configurations for (N). The topological susceptibility has smaller autocorrelation times and a 
correspondingly larger sample size. The better agreement for (Q 2 ) can be used to argue that 
the differences between the bonding boxes vanish with increasing sample size, as expected. 
The bottom table gives the dimensions, defined in Fig. 1111 for the different bonding boxes Bi. 
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of the quark overlaps becomes unstable because the overlap matrix becomes 
nearly degenerate. To render the evaluation of the fermionic interaction stable 
we need to decompose the overlap matrix into smaller non-degenerate blocks. 

To achieve this, we truncate the interactions to zero above a certain cutoff; it 
is a function of the pair sizes and determines which instantons interact with each 
other. Given this cutoff, we can built clusters of pairwise interacting instantons. 
The overlap matrix Tja is then decomposed into a direct sum according to these 
clusters. Not only does this render the numerical manipulations stable, but it 
also leads to a dramatic speed increase. Incidentally, we also use a cutoff for 
the gluonic interactions to get a better scaling with the number of instantons. 



4. Quenched simulation 

The main purpose of this section is to investigate how well the phenomeno- 
logical screening factor (|35[) can model the confinement/deconfincmcnt phase 
transition. To this end, we compare the topological susceptibility, for various 
values of the free parameters, to lattice data Normalising to the T = 
result, we get A = 206(8)MeV. 

We performed many simulations with different sets {T*. AT}. We report a 
few of these in Fig. 1171 which is representative of the general findings. Namely, 
that the screening factors do not capture the lattice data well. In particular, 
we found that the IILM result for the topological susceptibility decays too fast: 
the power-law-like behaviour, e.g. \ cx T~ s from the dilute gas approximation, 
at temperatures above T* is not compatible with the available lattice data. 
The screening factor d\- generically underestimates the lattice result for rather 
moderate temperatures. Only higher T* could remedy this; however, for such 
values the topological susceptibility is too high in the low temperature regime. 
This can also be seen in Fig. [17] for the second screening factor Sp\ it has a 
gentler switch-on of the plasma screening effects, and also leads to rather large 
values for the topological susceptibility at low temperature. 

We actually found that without any plasma screening, and therefore also 
in the low temperature region, the IILM produces a susceptibility that rises 
with temperature. This is not what the lattice predicts, namely an almost 
constant \- The quenched IILM at finite temperature does not seem to give a 
good description of the pure Yang-Mills sector: that it fails around the phase 
transition could be expected since the crude phenomenological screening terms 
can only be expected to be a rough effective description; however, we find that 
the IILM fails to qualitatively describe the pure gauge dynamics at low and high 
temperature. 



Given the fact that the dilute gas of non-trivial holonomy calorons 11] gave 
rather encouraging results, it might be that these are essential to model QCD 
in the quenched case. In particular, it will be interesting to find out whether 
these degrees of freedom can capture the smaller decay at high temperature or 
whether other degrees of freedom are dominant. 
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Figure 17: A striking feature of the quenched simulation is that the topological susceptibility 
decays generically too fast, as compared to lattice data. At low temperatures the IILM predicts 
a rising topological susceptibility, in contrast to the lattice data, which is roughly constant. 
The screening factor underestimates the lattice result, whereas d^, tends to overestimate 
it at temperatures around T„. Once the plasma screening is fully effective, both lead to too 
fast decays: the dilute gas approximation gives an approximate power-law of \ <x T~ 8 which 
is not compatible with the available lattice data. 
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5. Chiral symmetry restoration 



Compared to the quenched sector, full QCD has the additional property of 
(softly broken) chiral symmetry. Based on it, powerful analytical approaches 
like chiral perturbation theory have been developed. A strong point of the IILM 
is its rather accurate description of the chiral properties of QCD. As compared 
to the quenched sector, this additional structure can be expected to improve 
the agreement between the IILM and the lattice at finite temperature. 

We will repeat the analysis of the previous section for the unquenched IILM. 
Instead of the topological susceptibility, which we want to compute after all, we 
will use the chiral susceptibility [3J to estimate the parameters T* and AT. 
Note that there is a controversy with regard to the lattice results for the critical 
temperature 0, EH Q with differences on the order of 20 MeV. 

The chiral susceptibility is defined by 

X - qq = ^lnZqco = J J (H>(x)W(v)) ~ (J W>{x)\ , (38) 

with Zqcb the QCD partition function. It is actually composed of a connected 
and a disconnected part, and the latter is most sensitive to the chiral transition 
[20I ] . In terms of propagators both parts are given by 

JttS a (x,x)\ ) a - (lj TrS A {x,x)) A \ , (39) 

X% q = -([ TrS A (x,y)S A (y,x)) A +xi, (40) 

and the expectation values are over the gluon fields. These formulas make it 
clear that the disconnected part is the main order parameter of the transition. 
Within the IILM, the propagator is the sum of a low and a high frequency 
part; the latter is usually just approximated by the free massive propagator. 
As such it gives an infinite contribution to the disconnected part and has to be 
subtracted at T — 0. A finite contribution remains at T 7^ 0; it can be derived 
from the free partition function [18j and turns out to be small compared to the 
low-frequency contribution, and we will ignore it. This is expected, as the chiral 
transition is driven by the dynamics of the interacting instanton ensemble, and 
so the contribution from the free propagator should be negligible. By the same 
reasoning, we will ignore the cross terms between the low frequency and the 
free propagator. It is then straightforward to compute the chiral susceptibility 
in the IILM. The propagator, truncated to the finite dimensional low frequency 
part, is given by 

S A (x, y) = (p + m)- o l = ]T — L— e„(aO£t (j,) , (41) 

Til %/\ r fi 

n 

where A are the eigenvalues of the matrix defined in (fT0|) ; note that we used the 
fact that the set {£n} of zero modes is assumed orthonormal. The susceptibilities 
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in the IILM are then given by 




(42) 
(43) 



For notational simplicity we have omitted the subscript A on the expectation 
value, but the average is again over the gauge fields. Remembering that the 
non-zero eigenvalues are paired due to chiral transformations, we find that 



X 



m * — ' + A 2 ' 



\v\ s -^2(m 2 - \ 2 < 



X qq V ^ (m 2 + A 2 ) 2 ' Xqg ' 1 ' 

where \v\ is the number of zero modes and is related to the topological charge 
Q through the index theorem, v = Q. 

Now, the chiral susceptibility in Q is isospin symmetric, whereas we have 
broken it explicitly, i.e. m u ^ m^. In order to obtain effectively isospin sym- 
metric results, we will compute 

x* = 2ftiF' (46) 

with ifi = ^(m u + nid) the mean quark mass and m u /md = const. This leads 
to 

^Xqq = \™\x c uu + \ m \x C M + m u m d ((uu ~ (uu))(dd - (dd))) . (47) 

It turns out that the chiral susceptibility, with all plasma screening effects re- 
moved, peaks around T c w 120 MeV, and decays very slowly afterwards. Screen- 
ing effects in the IILM tend to increase the dilution of the ensemble, and thus 
they will be responsible for the faster decay at higher temperature. Inciden- 
tally, with screening effects included down to very low temperatures, below T c , 
the peak in the chiral susceptibility disappears completely; this gives a first 
constraint on T». 

The lattice data shows a peak at T c w 16OMe"V0 and no tuning of the 
parameters {T*, AT} will shift the chiral phase transition in the IILM towards 
the lattice result. 



5 The newer analysis in Q| results in a shift towards lower temperatures by about 5 MeV. 
This is not a significant change with regard to the IILM prediction of T c £3 120 MeV and we 
will ignore it in what follows. 
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Figure 18: The chiral susceptibility in the IILM predicts a chiral phase transition at 
T c w 120 MeV. Given the constraint that T* > 120 MeV this is a robust result. The dif- 
ferent parameter sets were chosen such that the plasma screening effects are effectively fully 
enhanced at about 160 MeV, the transition predicted by the lattice result (note that we use 
the continuum extrapolation of Q). The first curve (solid circles) falls slightly short of this 
expectation. Taking the errors at face value (including those for the lattice result) the IILM 
is off by roughly 3cr. Considering the newer analysis |2| that shifts the lattice result to lower 
temperatures, the agreement improves slightly. 



To fix the free parameters, we simply demand that the IILM result for the 
chiral susceptibility does not exceed the lattice result at high temperature, i.e. 
that the rapid decay due to screening effects should be around T « 160 MeV, and 
that its peak remains intact. We've computed the chiral susceptibility for four 
sets of parameters to estimate the systematics, see Fig. [T5] The large errors are 
due to the large systematics in A that is used to set dimensions. By construction, 
the data below T c ps 120 MeV is very robust as the screening effects are almost 
completely negligible, and the location of the peak remains fairly unaffected. 
Above T c the results do not depend too strongly on the parameters {T*, AT}, 
apart for the first curve (solid circles) which does seem to have too big a T*-value. 
Given the phenomenological terms, see ([35]). and the constraint T* > 120 MeV, 
the IILM predicts a chiral phase transition at T c s» 120 MeV. 

One of the successful predictions of the IILM has been chiral symm etry 



restoration based on the formation of instanton-anti-instanton molecules 15|, 
(l6| . [35| : through pairing up of instantons and anti-instantons, and the strongly 
localised quark wavef unctions, the Dirac operator spectrum develops a gap and 
the Casher-Banks relation tells us that the quark condensate vanishes, i.e. chiral 
symmetry is restored. 

Naively we would expect that the pair concentration is highest around the 
phase transition. This is, however, not what we find and, interestingly, the 
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Figure 19: Contrary to expectation, the pair concentration does not peak at the chiral phase 
transition, although it is large with 50%, but, interestingly dips at the phase transition. Also, 
it keeps growing to 70% at roughly T 300 McV; beyond that point it decays. 
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Figure 20: The instanton density compared to the dilute gas result. The data clearly shows 
that even for temperatures as high as T = 400 MeV the ensemble is not compatible with a 
dilute gas. The higher density in the IILM is due to the energy stored in pairs. 



number of pairs seems to drop at the phase transition. In fact, the concentration 
of molecules keeps growing beyond T c , see Fig. [19j and at T c only half of the 
instantons are paired up into molecules. The maximum concentration of 70% is 
reached at T w 300 MeV, beyond which it starts to decay rather slowly. Note 
that even at T = 400 MeV the system is still far from a dilute random gas of 
individual instantons, see Fig. 1201 It is, however, important to note that the 
identification of pairs is correlated to our specific definition of a bonding box, 
and even though the Monte Carlo results did not depend on the definition of 
the box, the number of instanton-anti- instantons does certainly depend on it. 

Despite the fact that the ensemble is distinct from a random gas, the quark 
condensate tracks the dilute gas approximation already for moderately high 
temperatures, see Fig. 1211 The reason is that the Dirac eigenvalues for pairs 
are large compared to the quark masses; this in turn follows from the fact that 
pairs have lined up along the time direction and thus their separation becomes 
smaller as the temperature increases. Therefore they are negligible compared to 
the zero modes. Bearing in mind the cluster decomposition of the quark overlap 
matrix Tja, the quark condensate is given by 



— £ <i« 

111 a L — ' 



clusters 



m a 



(48) 



Note that this has almost the form of the dilute gas result, for which Tja = 
and hence A n = 0, 

<99>Idga = ^T^- (49) 



in,. 
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Figure 21: Although at T 300 MeV the pair concentration is high, the quark condensate 
already tracks the dilute gas approximation (solid line). 



It turns out that (l^effl) has a temperature dependence that follows the dilute 
gas approximation for (N), see Fig. [22j The straightforward explanation that 
\v c s I is approximately equal to the number of unbonded instantons did not hold 
up to scrutiny, especially (^unbonded) does not behave as (N/V)\ DGA . Thus 
the population of unbonded instantons is still too large. But the clustering is 
defined through cutoffs in the interactions; these cutoffs can therefore also be 
used to identify a truly non-interacting sub-ensemble in the IILM. By definition, 
this population behaves according to the dilute gas approximation. We find that 
its density becomes important at T w 300 MeV which explains the behaviour of 
the quark condensate. 

We can conclude that, quite generally, the quark condensate is rather insen- 
sitive to the details of the molecule population. Presumably, only its concentra- 
tion determines the effective number of zero modes v e g, respectively the density 
of the non-interacting instantons. 

6. Topological susceptibility and axion mass 

6.1. Topological susceptibility 

The topological susceptibility approaches the dilute gas approximation for 
temperatures above approximately 250 MeV, and, like the quark condensate, 
does not seem to be sensitive to the high concentration of instanton-anti- 
instanton molecules, see Fig. [23] 

Guided by the ideas of [16[ that the IILM leads to a mixture of a highly 
correlated and a random component, it is natural to identify the unbonded 
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Figure 22: It turns out that the effective number of zero modes, |^ e ff|, follows the temperature 
dependence of the dilute gas approximation for the density (solid line). Together with the 
observation that the condensate is dominated by v a g at high temperatures, see text, this 
explains why the condensate has a temperature dependence that follows closely the dilute gas 
result. At temperatures around T m 300 MeV we can attribute this effect to a sub-population 
of non-interacting instantons. 
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Figure 23: The IILM reproduces the topological susceptibility in the dilute gas approximation 
for temperatures where the ensemble is still far from dilute. We found that this phenomenon 
cannot be attributed to a population of unbonded instantons. Rather, it can be explained 
by the population of non-interacting instantons, Fig. I22I which is smaller than that of the 
unbonded instantons. At higher temperatures the distinction between unbonded and non- 
interacting instantons becomes ever more unimportant. 

instantons as the dominant contribution to the topological susceptibility; the 
molecules, at least at zeroth order, do not lead to charge fluctuations and thus 
cannot account for the topological susceptibility. For a random ensemble N r it 
follows that 

X oc (N?) - (N r ) 2 oc (JV r ) . (50) 

We found, however, that the population of unbonded instantons is far too large 
to be responsible for the small topological susceptibility, under the assumption 
that N u b is a completely random ensemble, see Fig. [Ml Again, as in the case 
of the quark condensate, it is the non-interacting instantons, not the unbonded 
instantons, that fit the bill. Note that for higher temperatures, T ~ 600 MeV, 
where the pair concentration is still quite high, the concentration of unbonded 
and non-interacting instantons becomes equal and there is no more distinction 
between the two. For lower temperatures a coupling between unbonded and 
bonded instantons still exists. 

Note that, in the chiral limit, the topological susceptibility is related to the 
quark condensate through chiral perturbation theory. It is therefore consistent 
that the condensate and the topological susceptibility behave rather similarly. 

Finally, it is worth noting that the topological susceptibility, and also the 
quark condensate, are not overly sensitive to the free parameters T* and AT, 
see Fig. [25] 
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Figure 24: For moderately high temperatures, the population of unbonded instantons can- 
not be responsible for the topological susceptibility (solid line); we have checked that 
(JV^ 6 ) — (N n b) 2 ~ (JVufc)- The unbonded instantons still interact too strongly with the highly 
correlated instanton-anti-instanton molecules. For higher temperatures, the unbonded in- 
stantons become indistinguishable from the non-interacting instantons that can account for 
the topological susceptibility. The latter instantons play the role of the random sub-ensemble 
along side the 'molecular' ensemble, following the ideas of [la ]. 
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Figure 25: Both the topological susceptibility and the density are fairly constant for low 
temperatures, except for a very gentle growth from T = towards a maximum at around 
T = 100 MeV. In contrast to the quenched case, the IILM is not as sensitive to the free 
parameters, {T*,AT}, that effectively describe the phase transition. The first curve (solid 
circles) could be argued to be unphysical as it delays the decay in a rather unnatural fashion. 
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6.2. Temperature- dependent axion mass 

In order to solve the strong CP problem, Peccei and Quinn introduced a new 



field into the Standard Model of particle physics [32j, |3l| . It was soon realised 
that this field gives rise to a new light particle, a pseudo-Goldstone boson, the 
49l 46 1 . At the time it seemed natural to tie the axion to the electro- 



axion 



weak scale but laboratory experiments have ruled out such an axion. Leaving 
the scale of the PQ field f a free, and large to evade the previously mentioned 
constraints, the so-called 'invisible' axions were born, e.g. 0113, HE 



Through their weak couplings to ordinary matter, invisible axions can play 
the role of a dark matter candidate. They have a rich phenomenology, and can 
be produced through a variety of production channels: the thermal scenario 
(similar to WIMP production), cosmic string decay (axions are the the Gold- 
stone boson of a spontaneously broken C/pq(1) symmetry) and the so-called 
misalignment mechanism. The latter is essentially the PQ mechanism: due to 
the anomalous J7pq(1) symmetry an axion mass term is generated through the 
coupling of the axion to the topological charge 



This term combines with the vacuum angle 6 and the axion field can be shifted 
to get rid of 9, i.e. the vacuum angle becomes a dynamical field, the axion angle 
6 a . This is the key insight because now we can evoke the principle of least action 
to argue that 9 a — > dynamically to solve the strong CP problem. 

Indeed, integrating out the gluons, we can determine the axion mass from 
the effective action, 

exp(-FFeffW) = J[dA] det(#> + M) cxp(-S 9 - S a - g ) , (52) 

according to 



mi 



d 2 V cS 



0=0 



. (53) 



and the axion eventually evolves to its minimum at 9 a = [43j. It also demon- 
strates that the mass for the QCD axion is set by the topological susceptibility. 

Of all the production channels, the misalignment mechanism is most sensitive 
to the axion mass. From above it is clear that the axion mass is inherently a 
non-perturbative problem. At high temperatures, the dilute gas approximation 
to the instanton ensemble can be used. However, at lower temperatures it breaks 
down. From our determination of the topological susceptibility we can for the 
first time give a well-motivated axion mass that covers all temperatures down 
to T = 0. 

In Fig. [26] we display the axion mass together with a fit and its error range 
mainly due to the error in A. The data suggests that the axion mass turns 
into the dilute gas approximation rather quickly; around T ~ 300 MeV the 
differences are negligible, and the fit takes this into account. We have seen in 
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the previous sections that the IILM predicts a phase transition that is slightly 
too low. Using the lattice data for the phase transition temperature T c , we 
also include a tentative fit to what the ultimate lattice data might look like, 
and again we impose the dilute gas limit at moderately high temperatures. 
It is worth noting that we currently do not have lattice data available that 
could corroborate such a result; remember that in the quenched case the lattice 
data did not behave according to the dilute gas result at high temperatures! 
Given the state of lattice calculations it should not take too much longer before 
a comprehensive lattice study with physical quark masses and across a wide 
range of temperatures will be available to give the exact axion mass. 

Within the IILM some progress can be made by including the non-trivial 
holonomy calorons; a dilute gas study in the quenched sector has given en- 
couraging results [ll| that these degrees of freedom might play a role in the 
confinement /deconfinement transition. It would be very interesting to investi- 
gate their role in the unquenched sector, where we might expect less dramatic 
qualitative changes for chiral quantities, such as the topological susceptibility, 
because chiral properties are reasonably well modeled by the IILM. However, we 
certainly expect a closer agreement with the lattice results if the QCD vacuum 
is indeed dominated by non-perturbative, topological fluctuations. 

A second improvement with regard to finite temperature effects is the impli- 
cation of quantum interactions on the plasma screening effects. In particular, 
we have found within a toy-model that the topological susceptibility can change 
qualitatively if the screening effects become subdominant to the quark zero 
mode interactions. This would favour a higher concentration of molecules and 
therefore reduce the non-interacting instanton density that sets the axion mass. 
Specifically, it would lead to a faster decay of the axion mass. 

The IILM result has a slight rise at low temperatures. We do not have enough 
low-T data to determine the shape of this rise; however, in the gauge sector we 
found that the topological susceptibility had a roughly linear dependence on 
T. We will therefore constrain the fit at low temperatures to be a first order 
polynomial in temperature. It is given by 

mlf a = 1.46 1(T 3 A 4 1 + °- 50T/A 748 , < T < 0.45 , (54) 
1 + (3.53T/A) 7 ' 48 V ; 

where A = 400 MeV, and the errors will be mainly due to the uncertainties in 
A. 

Inspired by the lattice result for the chiral phase transition [i3|, we include 
a fit that delays the decay of the topological susceptibility until T c « 160 MeV. 
As mentioned above, there is a rather large disagreement between different 
lattice collaborations, and the chiral phase transition could occur at higher 
temperatures still. We assume again that the dilute gas limit will be recovered 
rather quickly above the phase transition. This assumption is hard to justify 
quantitatively given the lack of lattice data in that regime; the dilute gas result 
is really the best estimate we have at the moment. For the lattice inspired 
fit to smoothly connect with the dilute gas limit at around T rj 200 MeV, we 
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Figure 26: The mass for the QCD axion follows from the topological susceptibility, m?/^ = 
X- The fit goes over to the dilute gas approximation for moderately high temperatures, 
in accordance to the IILM data. Note that the large errors are mostly due to the large 
uncertainties in the determination of A, used to set dimensions. We also include a similar fit 
that is slightly shifted towards higher temperatures to mimic the phase transition as seen on 
the lattice, and a simple power-law approximation to the dilute gas limit (DGA), see Q59|l . 
We see that such a simple power-law approximation to the full result certainly has its merit 
as it is fairly accurate given the analytic simplicity. 
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patch together two different rational functions with the help of a washed-out 
step function, a(T) = | (l — tan T ^f 07 °' 40 ) . The result is 



mlfl = 1.46 1(T 3 A 



3 A 4 



1 + 0.3T/A 
1 + (2.5T/A) 8 ' 3 



o(T) 



1 + (3.4 T/A) 



7.4 



a(T)-l 



(55) 



for the temperature range < T < 0.41 GeV. 

These fits cover the low temperature regime, bounded by what we call Tdga- 
Given our aim to improve on the current axion mass computations, we will also 
be more systematic for high temperatures and take the effects of quark thresh- 
olds into account. To that end, remember that we really are using the language 
of effective field theory: the results are given as a function of g( 3 \ the strong 
coupling for Nf — 3 active flavours. It is well known that in order for S-matrix 
elements to be smooth, the 'free' parameters, such as masses and coupling con- 
stants, become 'discontinuous' [47, II, 3(1 33|; actually, the parameters are not 
discontinuous but belong to different theories that are matched at quark mass 
thresholds, i.e. g^ = g^ + 0((</ 4 ') 2 ). Since part of the quantum effects in 
the instanton density are two-loop improved, we expect that such 'discontinu- 
ities' do arise. In the present case the coupling in the exponential is the only 
place that displays two-loop improvements, and, at the threshold, the dilute gas 
approximations are related to each other by 



X (4) =X (3) exp( 7M ). 



(56) 



The factor 73 4 is essentially given by the Dirac determinant of the newly ac- 
tive quark flavour. In practice we just determine 73,4 such that the dilute gas 
approximations are smooth across the quark threshold. 

At high temperatures the dilute gas result is almost given by a power-law. 
We found that the corrections could be straightforwardly taken into account by 
a higher order polynomial in log-log space; therefore the fit takes the form of an 
exponential of a polynomial in logarithms. We find 
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with I — In -? , and the different parameters are given by 



N f 




4^ 


4^ 


4 Nf) 


3 


-15.6 


-6.68 


-0.947 


0.555 


4 


15.4 


-7.04 


-0.139 




5 


-14.8 


-7.47 


-0.0757 





(58) 



The jNr-i,Nf -factors have already been absorbed into d Q 
0.444 and 74,5 = 1.54. 



they are 73,, 
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We conclude by giving a very simple approximation to the dilute gas result 
in the form of a power-law, as in earlier work 41. 3|. 



2 «aA 4 

m « = imw> (59) 

where n — 6.68 and a — 1.68 10~ 7 , from (|55|) . and compare well with the more 
recent study [4J; it compares well with [4j. We believe it is a rather lucky 
coincidence that such a simple fit, based solely on the high temperature regime, 
still gives such a good overall approximation to the much more elaborate result 
of the IILM simulations, see Fig. [26l The small qualitative differences are that 
the power-law approximation: overshoots the IILM result at high temperatures, 
due to the wrong running of the QCD /3-function, and underestimates the axion 
mass at low temperature where it is by construction constant whereas it reaches 
the T = limit from above in the IILM. Given its analytic simplicity and its 
unexpectedly good agreement, such a simple power-law has certainly its merit. 
Whether such a conclusion pertains to an improved IILM including the more 
general calorons, and ultimately to the lattice, remains an open question. 



7. Conclusion 

We have been able to improve on the finite temperature interactions in the 
IILM. The numerical framework we set up in the first paper [3] could success- 
fully be implemented at finite temperature as well, and well-defined interactions 
that lead to a consistent thermodynamic limit have been derived. 

Using these improved interactions, we investigated the IILM at finite tem- 
perature with light, physical quark masses. For these small quark masses we 
have found that the usual 'random' Monte Carlo sampling is very inefficient and 
can even break down. Using the results form [45j |. where we introduced biased 
Monte Carlo techniques and, in particular, adapted the Unbonding-Bonding 
algorithm to the grand canonical ensemble, including a biasing scheme to deal 
with the orientation-dependent interaction of the IILM, we could run efficient 
simulations at finite temperature. We have found that the screening factors, 
from single instanton quantum fluctuations, will lead to a dilute, random en- 
semble at high enough temperatures. There is, however, a possibility that this 
trend could be reversed when more quarks become active; this does seem rather 
unnatural though. We want to point out that the high-frequency quantum inter- 
actions from overlapping instantons might be important to settle this question. 

In the Yang-Mills sector we found rather poor overlap of the IILM and the 
lattice data. Most strikingly, we found that the IILM failed to reproduce the lat- 
tice data, even in a qualitative manner, in the regions where the phenomenologi- 
cal factors that mimic the phase transition in the IILM are unimportant, namely 
at low and high temperatures. To model the pure gauge sector, the IILM might 
have to be generalised to include the non-trivial holonomy calorons. Ultimately 
this is needed in any case because, as for the zero temperature IILM, confinc- 
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ment is still lac king and these more general degrees of freedom, the KvBLL 
calorons 2^, 2|| 24 , 2|| , might play a cruciaH. 

In the unquenched sector, we investigated the chiral susceptibility and the 
quark condensate to gauge the free parameters introduced by the phenomeno- 
logical screening factors; we found that the IILM is not overly sensitive to these, 
given some mild restrictions. The IILM unambiguously leads to chiral symmetry 
restoration at T c = 120 MeV, slightly too low as compared to lattice data. 

We investigated the population of instanton-anti-instanton molecules and 
found that, rather surprisingly, the maximum concentration does not occur at 
T c but rather at higher temperatures. A large population of instanton pairs 
prevailed to fairly high temperatures, indicating that the dilute gas limit of the 
instanton ensemble is only reached far beyond the critical temperature. How- 
ever, the quark condensate and the topological susceptibility behave according 
to the dilute gas result much earlier. We could attribute this behaviour to a 
sub-ensemble of non-interacting instantons, distinct from the unbonded instan- 
tons that still interact considerably with the instanton-anti-instanton molecule 
population. At higher temperature, the distinction between unbonded and non- 
interacting instantons becomes irrelevant. 

Given the topological susceptibility, we have presented a fit to the axion 
mass. We paid due attention to extrapolate the axion mass to higher tempera- 
ture by including threshold effects due to heavier quarks that become active as 
the temperature rises. The main improvement, however, is a real computation of 
the low temperature axion mass that matches smoothly to high temperatures. 
Considerations of the anthropic axion for which 9 and f a are considered free 
parameters, one 'environmental' and the other fundamental, need knowledge 
of the axion mass for all temperatures, which this work provides. Comparison 
with lattice data leaves open the intriguing possibility that the high temper- 
ature axion mass does not behave according to the dilute gas result based on 
Harrington-Shepard calorons. A considerably different fall-off behaviour of the 
axion mass will change the cosmological bounds decidedly: the classic axion 
misalignment scenario, where 9 a is set by its rms fluctuations at the time of 
symmetry breaking, would get a weaker upper bound. This conclusion relies on 
our findings from the pure gauge simulations, where the lattice topological sus- 
ceptibility did not fall off as quickly as predicted by the IILM. The unquenched 
case is, however, significantly different and chiral symmetry, which has proved 
powerful at zero temperature, might constrain the discrepancies between the 
lattice and the IILM. Especially, a higher instanton-anti-instanton molecule 
density, possibly due to the non-trivial holonomy calorons or weaker screening 
effects for the strongly overlapping pairs (the effective size of such a pair will 
be smaller), would lead to stronger upper bounds and could possibly rule out 
the classic axion window. This was one motivation for the present work, but 
within the present IILM this expectation could not be corroborated. Further 



"Recent lattice studies see evidence of the lumpy structure characteristic for an ensemble 
of these new caloron solutions |17| . 
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investigations are necessary to settle this issue. It might turn out, however, that 
the non-interacting sub-population that sets the axion mass in the IILM will 
not change considerably; such robustness within the IILM could be interpreted 
as evidence that the axion mass is rather insensitive to the details of the QCD 
phase transition. Given the advance of lattice QCD simulations, the ultimate 
axion mass determination might be available in the near future. 
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Appendix A. Gluonic Interactions 

The ratio ansatz for an instanton-anti-instanton pair is defined by 



where O — 0*C>2, with 0% the respective colour embeddings. A global colour 
rotation has been performed to bring the gauge potential into this form. Since 
we are ultimately interested in the action, we do not need to bother about it 
since the action is gauge invariant. Instanton-instanton and anti- instanton- 
anti-instanton pairs differ by having either only fj or r\ in the above formula. A 
brute force computation then gives 



K»Kv = I + ( Tr0t ° + (vOv)^)J + (fjOv) P ^i^ 

+ {r]O t Ovi) ll py a J^p V( y + (r)Or)) a p, ap {r]Or])rjvj3< T Kpp l/l j ■ (A. 3) 




TfcAIIiOc, {an, pi}) + O a % v dMx, {x 2 ,p 2 }) 



,(A.l) 



n(x,{y,p}) 



1 + IIi (or, {xi,pi}) +n 2 (x, {x 2 ,/9 2 }) 
np 2 shrill 

pr cosh 2£I - cos ^ ' 



(A.2) 
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The different terms have the following form 

i = (1 + ni 4 +n2)2 [{d^u^d^Ui) + (d„d v ii 2 )(d^d v u 2 )} 

• [(^^n 1 )(^n 1 )(^n 2 ) + (^a„n 2 )(9 Al n 2 )(^n 1 ) 



(i + n! + n 2 )3 

+2{a lt d v n 1 )(d li n 1 ){d v n 1 ) + 2(^^n 2 )(^n 2 )(^n 2 )] 

+ (1 + ni 4 +n2)4 ^n^no^n^n,) + 3(a At n 2 a M n 2 )(a A1 n 2 a Al n 1 ) 

+3(9 Al n 1 a M n 1 ) 2 + 3(^n 2 ^n 2 ) 2 + 2(^n 1 ^n 1 )(^n 2 d / ,n 2 ) 

+(a Al n 1 a AJ n 2 ) 2 ] . (a.4) 



3 = (i + n 1 + n 2 )4 (^ni^ni)(^n2^n 2 ) • (A.5) 



^ = (i + n 1 + n 2)2 ( ^ ni)( ^ n2) 

+ (1 + ni 4 + n2)3 [(d^n^id^dM + (a Ai a„n 2 )(a ff n 1 a a n 1 ) 

-2(a M a (T n 1 )(a,n 2 )(a (T n 2 ) - 2(d li u 1 )(d tr u 1 ){a v d a u 2 ) 
-2(^5 CT n 1 )(a CT n 1 )(9 I/ n 2 ) - 2(^n 1 )(d 1 ^n 2 )(d CT n 2 )] 

+ ( i + ni 4 +n 2 )4 [-(^no^no^n,^^) 

-(^n 2 )(a„n 2 )(a CT n 1 9 CT n 1 ) + s^n^^n^n^n!) 
+3(^n 1 )(^n 2 )(a ff n 2 a CT n 2 ) + 3(^n 1 )(a y n 2 )(9 ff n 1 a CT n 2 )] . (a.6) 



4 

(i + n7+ n 2 ) 2 



+ (1 + ni + n2)3 [(d^Kdpd^KdM + (^n 1 )(a,9,n 2 )(a (T n 1 )] 

+ (i + ni 8 + n 2 )4 (Wid^id^id^) . (A.7) 



V- = (1 + ni + n2)4 (^n 1 )(9 p n 2 )(^n 1 )(9 ff n 2 ) . (A.8) 



K » pv ° = (i + ih + n 2 )4 (^n 1 )(^n 2 )^n 1 )(a g n 2 ) . (A.9) 
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Appendix A.l. Exact Interactions 

When computing the look-up tables we use global translations and rotations 
in R 3 , and periodicity in S 1 , to place one instanton at the origin and the partner 
at y' 3 = R s = y/R~R~i = lyl 1 - yl 2 \ and y' 4 = R t = |y t 71 - y[ 2 \- f are the 
instanton centres. The rotation will reemerge in contractions of Ri with the 
colour structure as we will now see. The relation between the spatial position 
vector Ri and i?- = (0,0, R s ) is given by the following 0(3) rotation matrix 

R'i = OljRj, (A.10) 

O l3 = |i, (A.11) 

and the other components of the rotation matrix are irrelevant. 

Note that with the choice of R[ the integrands are 0(2) symmetric in the 
subspace orthogonal to the 3-direction. Denoting the spatial arguments of the 
't Hooft potentials, Tl(x, {y, p}), by Xi and Xi = Xj — Ri, we can extract the 
Ri from the integrands with help of the following formulas, which we order 
according to the tensor structure in the a^-dependence of the integrand 

J x, = O a J x' 3 . (A.12) 
dij f x? + O l3 O j3 f (x' 3 2 - xf) . (A.13) 



J Xii 

J XiXjXk = (S tJ O k3 + 6 kl O : j 3 + S lk O t3 ) J xfx' 3 (A.14) 

+ O l3 O j3 O k3 J (x'i - 3xi 2 x^) . (A.15) 

J XiXjX k x h = {Sijdhh + SihSjh + SihShj) J xfx'2 (A.16) 



(A.17) 



+ {5 l0 O k3 O h3 + perm.) J {x'fx? - xfx' 2 2 ) 

+ O l3 O j3 O k3 O h3 J (x'i - 6xfx'i + Sx'fx'i) . (A.18) 

Insertion of x can be constructed from these. Incidentally splitting the dif- 
ferent integrands according to the above formulas is the most stable procedure 
numerically. Taking into account the antisymmetry of the 't Hooft symbols, we 
end up with the following integrands 
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(A.19) 



Note that to achieve good numerical precision, we need to subtract the one- 
instanton integrands from the above before performing the numerical integration 



(i + n! + n 2 ) 4 



((n' 1 ) 2 + (n 1 ) 2 )((n 2 ) 2 + (n 2 ) 2 ). 



(A.20) 
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Appendix A. 2. Asymptotic Interactions 

As explained in the main text, the small separation asymptotic formulas 
get contributions which have the same functional form as those for the large 
separation asymptotics but with different explicit integration limits, see Fig. [§1 
We will therefore start with the large separation formulas and not evaluate the 
integrals explicitly. 

Appendix A. 2.1. Large Separation 

As explained in Fig. [TJ the 't Hooft potential is approximated by the T = 
BPST form for small arguments up to = ^jx^x^ = /3/2; for larger arguments 
we use the r 3c i —¥ oo asymptotic form of the T ^ potential. These simple 
rational expressions allow us to factor out completely the 't Hooft potential of 
the partner instanton, whose argument is kept fixed, except for its contribution 
to the integration limit. For fixed J 2 , and depending on the integration region, 
this leads to the two integration limits 

9 1 + n 2 , 

Pi 

i + n 2 

z 3 d = — ; r 3d . (A.60) 

Pi 

We have defined p = np 2 j ' (3. In the following we will encounter the integrals 

r = io &mi L=JT- 

Finally, for r 3( i > ft/2 the approximate 't Hooft potential is independent of 
t. This will lead to some indices only running over the spatial set i = {1,2,3} 
whenever we use roman letters. 

The integral over / contains terms that do not mix the 't Hooft potential Hi 
and n 2 except for the denominators. Just as in the T — case, these can be 
transformed to exactly match the single instanton contributions by exploiting 
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scale invariance and independence of the action on /?. Remembering that we 
actually subtract the one-instanton contributions to get the interactions, we can 
neglect these terms altogether. We then end up with the following formulas 



fr 7 y; W! p d sds 

+ ( 32 ^ (fr# + i^tiw ) L (sTiy + sym ■ (A - 61) 

r 2 2 ^n^n r z ^ s^ds 
J J = W7rp aTWJ WTW 

+ 8 ^(ft^I M (7TTF + sym - (A - 62) 

/"r Ifi 2 2 g^jl f z ^ s^ds 

i v = 16l "(iWl FTIF 

- r P ""(1 + 11)3 + 8 - p (i + n) 3 J y J^TW 
+ l 16 ^(TTnF " I^wj L (s + W 

- " Pf + ^ (1 + 11)3 

(g^n)(gji) \ /• s 2 <fe 
-^TTW J 1, (sTTF + sym • (A - 63) 

At zeroth order, partial integration and the antisymmetry of the 't Hooft 
symbols can be used to simplify 

w ~* (i + ni 8 + n 2 )4 (^ni)(gpn 2 )(ft,ni)(a g n 2 ) . (A.64) 

With asymptotic behaviour 

y vp- io7r p v (i+n)3 y ( S 2 + i)4 

32 „ r (<9,,n)(<9 CT n) /■ s 2 rf s 

+ Y^irhrL —* +sym - (A - 65) 

/ t A^nH (W(^n) 

y w ^ ~ p ^ (i + n)3 y ( s 2 + 

8 (dpTlYd^m f s 2 ds 
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For K ppua no 't Hooft symbols can be used to exchange the index pairs 
(/i, v) f-> (p, a), and so we cannot simplify with a symmetry argument anymore 



f „ _ , 2 2 c (5 P n 2 )(^n 2 ) /•«? 
7 ww ~ Pl<v (i + n 2 )3 y ( s 2 + i 
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s 2 ds 
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Appendix A. 2. 2. Small Separation 

As explained in Fig. [51 the small separation asymptotic formulas get con- 
tributions from the large asymptotics. Also, in this case we have performed a 
global translation so that the instantons sit at zLR^/2. For the terms in J the 
integration limit is given by z 2 = ^P-{R/2) 2 . The T ^ terms given by the 
J integrals are to be interpreted as 8h(R — (3/2) J , i.e. they only contribute if 
the separation R is bigger than j3/2. In practice, these terms do not contribute 
because they are covered by the look-up tables. The proper small separation 
asymptotic formulas, that encode the repulsion through the gauge singularity, 
arc then given by the T = formulas, which we repeat here for convenience. 

Introducing another explicit upper limit z 
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and approximating the arguments x p ± R^/2 



we arrive at 
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Appendix B. Fermionic Interactions 

The Dirac overlap matrix elements are given by 

Some straightforward algebra leads to 
j = ^ 

(i + n / + n A )(i + n / ) 1 /2(i + n A )i/2 



i + n 7 v p ^ i + nj y V i + n,4 



i + n A /V v i + n 7 



J/3Aia (i + n / + n A )(i + n / )i/2(i + n A )i/2 
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(i + n / + n j4 )(i + n / ) 1 /2(i + n A ) 1 /2 

(d, X Ad a Il A ) {d, Xl - &5^) . (B.4) 

Appendix B.l. Exact Interactions 

The three contributions give rise to two different colour structures which can 
be combined with the help of the colour four- vector iup = iTr({7rj~), used for 
instance in 
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Appendix B.2. Asymptotic Interactions 
Appendix B.2.1. Large Separation 

At T ^ the quark zero mode has an additional factor \ °t exp(— r//3). 
Therefore we completely neglect the contributions from r > (3/2, and recover 
the T — formulas. For convenience we display the final results again here; 
details are given in The only difference with the T = formulas is that 

the upper limit z oc /3/2 instead of infinity 



2 2 IU^Ha 
8tt pj 



2 



{l+n A )V*J («2 + i)7/: 

, , ggn, p s 4 rfa 

" 87F ^(1 + ^)3/2 7 ( S 2 + 1)5/2- I 6 - 8 ) 

At zeroth order, we have that J Jp^ a — J Kp^ a = 0. For Jf)^ a this follows 
from the fact that within the r < (5/2 region we approximate \ ~* n and so 
the integration over Ij yields zero because of the anti-symmetry of the 't Hooft 
symbols. Integration over I a vanishes too because of 0(4) symmetry. 

Appendix B.2.2. Small Separation 

At zeroth order, i.e. x^±R^/2 — > x^, the contribution to Ip vanishes because 
of 0(4) symmetry. It turns out the large separation asymptotics falls off too 
quickly as R — > 0. However this is not important because in this regime the 
gluonic interaction is dominant. 
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